//Main.java
/**
 Main function of the Monte Carlo Ising simulation 
 **/

package simulation;

import static common.FSave.fSave;
import static common.PathJoin.pathJoin;
import static common.Print.println;
import static common.Range.range;
import static common.Todouble.todouble;

import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Date;

import org.apache.commons.math.stat.StatUtils;
import org.jdom.Document;
import org.jdom.Element;
import org.jdom.JDOMException;
import org.jdom.input.SAXBuilder;

import common.Stat;
import common.String2Array;

public class Main {
	static double Bi = 0;
	static double Bf = 0;
	static double Bs = 0;
	static double Ti = 0;
	static double Tf = 0;
	static double Ts = 0;
	static double dT = 0;
	static double J = 0;
	static int N = 0;
	static int MCS = 0;
	static int ThermalMCS = 0;
	static int[][] neighbors;
	// static String modelName = "two_d_square";
	static String modelName = "two_d_triang";

	// static String modelName = "one_d_chain";

	public static void main(String[] args) throws IOException {
		println();

		// read modelname
		// modelName = Cwd.argument(args, "-m");

		// read parameters from configuration file with JDOM
		SAXBuilder builder = new SAXBuilder();
		File xmlFile = new File("src/config/main.xml");
		Document document = null;

		try {
			document = (Document) builder.build(xmlFile);

		} catch (IOException io) {
			System.out.println(io.getMessage());
		} catch (JDOMException jdomex) {
			System.out.println(jdomex.getMessage());
		}

		Element rootNode = document.getRootElement();
		// reading general parameters
		Element general = rootNode.getChild("general");
		final double J = Double.parseDouble((general.getChildText("J")));// coupling
		// constant
		String method = "Glauber";
		if (J > 0) {
			method = "Wolf";

		}

		// magnetic field range
		final double Bi = Double.parseDouble((general.getChildText("Bi")));//
		final double Bf = Double.parseDouble((general.getChildText("Bf")));
		final double Bs = Double.parseDouble((general.getChildText("Bs")));

		if (Bi > Bf) {
			println("Final Magnetic field should be less then initial");
			System.exit(1);
		}

		// Temperature range
		final double Ti = Double.parseDouble((general.getChildText("Ti")));
		final double Tf = Double.parseDouble((general.getChildText("Tf")));
		final double Ts = Double.parseDouble((general.getChildText("Ts")));

		if (Ti > Tf) {
			println("Final temperature should be less then initial");
			System.exit(1);
		}

		// Delta tau
		final double dT = Double.parseDouble((general.getChildText("dT")));

		// Number of sites in spatial direction
		final int N = Integer.parseInt((general.getChildText("N")));

		/*
		 * Variables MCS and ThermalMCS should be calculated with respect to
		 * parameters B, N, J, dT, T
		 */
		// Number of Sweeps
		// final int MCS = Integer.parseInt((general.getChildText("MCS")));

		// Number of Sweeps
		// final int ThermalMCS =
		// Integer.parseInt((general.getChildText("ThermalMCS")));

		// neighbors array describing unit cell of the lattice. For any
		// model it will be 2d array
		Element model = rootNode.getChild("model");

		String tneighbors = model.getChild(modelName).getChildText("neighbors");
		neighbors = String2Array.Str2int2(tneighbors);

		// number of spatial dimensions
		int sp_dim = neighbors[0].length;

		int maxL = (int) (200 / dT);

		double[] tRange = new double[maxL - 1];

		for (int L = 2; L <= maxL; L++) {
			tRange[L - 2] = 1 / (dT * L);
		}

		// read parameters from command line

		// start loop for magnetic field
		for (double B : range(Bi, Bf + Bs, Bs)) {
			// start loop for temperatures
			double[] ttrange = range(Ti, Tf + Ts, Ts);
			if (B != 0) {
				ttrange = tRange;
				// ttrange = new double[1];
				// ttrange[0] = 0.5;
			}
			for (double T : ttrange) {
				java.util.Date Idate = new Date();
				// create lists in which we'll save data
				ArrayList<Double> tM = new ArrayList<Double>();// data for
																// magnetic
																// moment
				ArrayList<Double> tM2 = new ArrayList<Double>();// data for
																// magnetic
				// moment

				ArrayList<Double> tM4 = new ArrayList<Double>();// data for
				// magnetic
				// moment

				ArrayList<Double> tE_J = new ArrayList<Double>();// data for
																	// Energy J
				ArrayList<Double> tE_B = new ArrayList<Double>();// data for
																	// Energy B

				ArrayList<Double> tE = new ArrayList<Double>();// data for
				// Energy J + Energy B

				// start loop for magnetic field

				/*
				 * In this part we try get variable sweepStep. If we wont get it
				 * from file we will generate test simulation, find correlation
				 * length and it will be our sweepStep
				 * 
				 * From sweepStep we can find MCS and Thermal MCS
				 */

				CorrLength cor = new CorrLength(N, J, B, T, dT, modelName,
						method);

				int sweepStep = cor.SweepStep();
				int MCS = cor.MCS();
				int ThermalMCS = cor.ThermalMCS();

				// We need find number of layers
				int L = 1;

				// switch Richard's system
				// Not all temperatures will be accepted.
				if (B != 0) {
					L = (int) Math.round(1 / (dT * T));
				}
				// L = (int) Math.max(1 / (dT * T) + 1, 2);
				// for (String arg : args) {
				// if (arg.equals("-info")) {
				println();
				println("J: " + J);
				println("B: " + B);
				println("T: " + T);
				println("L: " + L);
				println("dT: " + dT);
				println("N: " + N);
				println("Correlation time: " + sweepStep);
				println("MCS: " + MCS);
				println("ThermalMCS: " + ThermalMCS);
				// }
				// }
				if (MCS == 0) {
					System.exit(1);
				}
				Lattice Lat = new Lattice(N, J, B, T, L, neighbors, method);
				// start loop for sweeps
				for (int sweep = 0; sweep < MCS; sweep++) {
					// update lattice
					Lat.update();
					if (sweep > ThermalMCS) {
						if (sweep % sweepStep == 0) {
							// record data at each appropriate sweep
							double E_J = Lat.Energy_J();
							double E_B = Lat.Energy_B();

							tE_J.add(E_J);
							tE_B.add(E_B);

							tE.add(E_J + E_B);
							double M = Lat.M() / (Math.pow(N, sp_dim) * L);
							tM.add(M);
							double M2 = Math.pow(M, 2);
							tM2.add(M2);
							double M4 = Math.pow(M, 4);
							tM4.add(M4);
							// double M6 = Math.pow(M, 6);
							// tM6.add(M6);
						}
					}
				}
				// save data to file

				fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
						+ N, "dT" + dT, "Energy_J"), "T" + T, (Stat.mean(tE_J
						.toArray(new Double[0])) / (Math.pow(N, sp_dim) * L)));

				fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
						+ N, "dT" + dT, "Energy_B"),
						"T" + T,
						Stat.mean(tE_B.toArray(new Double[0]))
								/ Math.pow(N, sp_dim) / L);

				fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
						+ N, "dT" + dT, "M"), "T" + T,
						StatUtils.mean(todouble(tM)));

				// fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
				// + N, "dT" + dT, "M6"), "T" + T,
				// StatUtils.mean(todouble(tM6)));

				fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
						+ N, "dT" + dT, "M4"), "T" + T,
						StatUtils.mean(todouble(tM4)));

				fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
						+ N, "dT" + dT, "M2"), "T" + T,
						StatUtils.mean(todouble(tM2)));

				fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
						+ N, "dT" + dT, "X"), "T" + T,
						StatUtils.variance(todouble(tM)) / T);

				fSave(pathJoin("results", modelName, "J_" + J, "B_" + B, "N"
						+ N, "dT" + dT, "C"), "T" + T,
						StatUtils.variance(todouble(tE)) / T / T);

				java.util.Date Fdate = new Date();
				// for (String arg : args) {
				// if (arg.equals("-info")) {
				println("Total waste of time - "
						+ (Fdate.getTime() - Idate.getTime()));
				// }
				// }

			}
		}

	}
}